# read all .csv of Locations and convert to spatial object
ff <- list.files("../Data/NZ_SRW_Argos_data/Download_13April22/All_NZ_locations/", full.names = T)
# convert to sf object
raw_tracks <- ldply(ff, read.csv) %>%
st_as_sf(coords = c("Longitude", "Latitude"), crs = 4326)
# reintegrate lat lon in columns because we will need them later
raw_tracks[c("Longitude", "Latitude")] <- st_coordinates(raw_tracks)
# project spatial points to an antarctic coordinate system
raw_tracks <- st_transform(raw_tracks, proj_antarctic) Remove satellite locations with quality Z and apply speed filter
# how many raw locations do we have
raw_tracks %>%
group_by(DeployID) %>%
dplyr::summarize(nb_locations = n())## Simple feature collection with 18 features and 2 fields
## Geometry type: MULTIPOINT
## Dimension: XY
## Bounding box: xmin: -1990937 ymin: -5204887 xmax: 6114844 ymax: 966054.4
## CRS: +proj=stere +lat_0=-90 +lat_ts=-71 +lon_0=0 +k=1 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
## # A tibble: 18 × 3
## DeployID nb_locations geometry
## <chr> <int> <MULTIPOINT [m]>
## 1 212499 1751 ((1042776 -4331716), (1042878 -4331340), (1042958 -…
## 2 212500 101 ((522744.7 -4335400), (527629.6 -4302724), (527647.…
## 3 215258 2131 ((1042341 -4321852), (1047942 -4337869), (1048065 -…
## 4 215259 1046 ((-191707.8 -4113366), (930421.1 -4271630), (930693…
## 5 215260 2 ((390322.3 -4915529), (2000106 -5204887))
## 6 215261 1170 ((1019345 -4226342), (1023114 -4316903), (1023437 -…
## 7 215262 522 ((687454.5 -4148217), (702342.2 -4128871), (705314.…
## 8 215263 2267 ((913058.8 -4310695), (919566.4 -4399070), (949488 …
## 9 46633 2670 ((934563.9 -4652730), (939939.8 -4656590), (940692.…
## 10 46635 1204 ((1034519 -4391005), (1035491 -4393757), (1036133 -…
## 11 46950 3983 ((-1422054 -4268807), (-1418681 -4263496), (-141107…
## 12 46955 1797 ((353001.1 -3093163), (354964.3 -3149689), (355070.…
## 13 Bill Wiremu 8359 ((-1990937 -2646971), (-1111275 -5027211), (814011.…
## 14 Rima 1850 ((1008044 -4357885), (1008238 -4366671), (1008510 -…
## 15 Rua 243 ((1012690 -4362642), (1012977 -4365951), (1017595 -…
## 16 Tahi 2969 ((-1892230 -5009731), (837729.4 -4381394), (890273.…
## 17 Toru 3560 ((1025301 -4338445), (1026972 -4339054), (1027273 -…
## 18 Whaa 1660 ((753046.5 -4163074), (761742.7 -3922445), (761757.…
# tag 215260 only has two locations so I am removing it
raw_tracks <- raw_tracks %>% filter(DeployID != "215260")
# filter z positions
raw_tracks <- raw_tracks %>% filter(Quality != 'Z')
# transform Date into POSIX format
raw_tracks <- raw_tracks %>%
mutate(Date = as.POSIXct(Date, format = "%H:%M:%S %d-%b-%Y"),
month = as.numeric(strftime(Date, format = "%m")))
# prefilter positions with the speed filter (see Reisinger et al. 2021 who did that before actually fitting the ssm)
# vmax in meters / s
raw_tracks <- ddply(raw_tracks, ~DeployID, function(d){
d$argosfilter <- sdafilter(lat = d$Latitude,
lon = d$Longitude,
lc = d$Quality,
dtime = d$Date, vmax = 5)
return(d)
})
raw_tracks <- raw_tracks %>%
filter(argosfilter != "removed") %>%
dplyr::select(-argosfilter)
# who is left? how many positions?
raw_tracks %>%
group_by(DeployID) %>%
dplyr::summarize(nb_locations = n())## # A tibble: 17 × 2
## DeployID nb_locations
## <chr> <int>
## 1 212499 1550
## 2 212500 95
## 3 215258 1829
## 4 215259 903
## 5 215261 1008
## 6 215262 430
## 7 215263 1980
## 8 46633 2413
## 9 46635 1053
## 10 46950 3486
## 11 46955 1679
## 12 Bill Wiremu 7330
## 13 Rima 1443
## 14 Rua 207
## 15 Tahi 2507
## 16 Toru 2989
## 17 Whaa 1472
g <- gbase +
geom_sf(data = raw_tracks, aes(col = DeployID, geometry = geometry), size = 0.8)
ggsave(g, file = "./Outputs/filtered_tracks_NZ.png", dpi = 300, width = 8, height = 8)
gFigure 1.1: Filtered satellite tracks of right whales tagged in the Auckland Islands in 2020 and 2021.
Assess mean time step between locations
timelaps_df <- ddply(raw_tracks, ~DeployID, function(d){
d$timelaps_hours <- NA
for (i in 2:nrow(d)){
d$timelaps_hours[i] = difftime(d$Date[i], d$Date[i-1], units = "hours")}
return(d)
})
# mean time step between locations (in hours)
a <- aggregate(timelaps_hours~DeployID, timelaps_df, mean)
a## DeployID timelaps_hours
## 1 212499 1.6367726
## 2 212500 4.8701537
## 3 215258 2.9402048
## 4 215259 2.8828732
## 5 215261 1.7521861
## 6 215262 15.0779429
## 7 215263 2.0706410
## 8 46633 1.5360574
## 9 46635 1.9188844
## 10 46950 1.8284794
## 11 46955 0.9882587
## 12 Bill Wiremu 1.2479754
## 13 Rima 1.2587930
## 14 Rua 4.6188660
## 15 Tahi 3.0898524
## 16 Toru 1.3542783
## 17 Whaa 0.9798064
mean(a$timelaps_hours)## [1] 2.944237
ggplot(timelaps_df, aes(timelaps_hours)) +
geom_histogram(binwidth = 5, col ="white", na.rm = T) +
theme_bw() +
facet_wrap(~DeployID, scales = "free") +
xlab("Time laps between successive locations (in hours)")Figure 1.2: Interval between positions recorded by the tags (in hours). These histograms are used to assess potential gaps in the transmission and to set the time step in the interpolated tracks.
source("Fun_TrackInterruption_SRW.R")
# need a time column to run
raw_tracks$time <- raw_tracks$Date
# run the custom function to split tarcks interrupted for more than 144 hours = 6 days
raw_tracks <- Fun_TrackInterruption(dataframe = raw_tracks, hours_gap = 144, id = "DeployID")
# set first timelaps of each segment to NA
raw_tracks <- ddply(raw_tracks, ~segmentid, function(d){
d$timelaps[1] <- NA
return(d)
})
# print number of locations per segments
raw_tracks %>%
group_by(segmentid) %>%
dplyr::summarize(nb_locations = n()) %>%
print(n = Inf)## # A tibble: 27 × 2
## segmentid nb_locations
## <chr> <int>
## 1 212499 1550
## 2 212500_a 84
## 3 212500_b 11
## 4 215258_a 1727
## 5 215258_b 102
## 6 215259_a 889
## 7 215259_b 14
## 8 215261_a 971
## 9 215261_b 37
## 10 215262_a 376
## 11 215262_b 30
## 12 215262_c 7
## 13 215262_d 17
## 14 215263 1980
## 15 46633 2413
## 16 46635 1053
## 17 46950 3486
## 18 46955 1679
## 19 Bill Wiremu 7330
## 20 Rima 1443
## 21 Rua_a 196
## 22 Rua_b 11
## 23 Tahi_a 2195
## 24 Tahi_b 2
## 25 Tahi_c 310
## 26 Toru 2989
## 27 Whaa 1472
# remove segments with less than 15 locations
toremove <- raw_tracks %>%
group_by(segmentid) %>%
dplyr::summarize(nb_locations = n()) %>%
filter(nb_locations <= 15)
raw_tracks <- raw_tracks %>% filter(!(segmentid %in% toremove$segmentid))
# print number of locations per segments
raw_tracks %>%
group_by(segmentid) %>%
dplyr::summarize(nb_locations = n()) %>%
print(n = Inf)## # A tibble: 22 × 2
## segmentid nb_locations
## <chr> <int>
## 1 212499 1550
## 2 212500_a 84
## 3 215258_a 1727
## 4 215258_b 102
## 5 215259_a 889
## 6 215261_a 971
## 7 215261_b 37
## 8 215262_a 376
## 9 215262_b 30
## 10 215262_d 17
## 11 215263 1980
## 12 46633 2413
## 13 46635 1053
## 14 46950 3486
## 15 46955 1679
## 16 Bill Wiremu 7330
## 17 Rima 1443
## 18 Rua_a 196
## 19 Tahi_a 2195
## 20 Tahi_c 310
## 21 Toru 2989
## 22 Whaa 1472
# now what's the average time between locations?
a <- aggregate(timelaps~segmentid, raw_tracks, mean)
a## segmentid timelaps
## 1 212499 1.6367726
## 2 212500_a 2.2604284
## 3 215258_a 2.2659622
## 4 215258_b 2.8483443
## 5 215259_a 2.6024343
## 6 215261_a 1.4457380
## 7 215261_b 5.1853704
## 8 215262_a 2.3169081
## 9 215262_b 5.2023563
## 10 215262_d 9.8943056
## 11 215263 2.0706410
## 12 46633 1.5360574
## 13 46635 1.9188844
## 14 46950 1.8284794
## 15 46955 0.9882587
## 16 Bill Wiremu 1.2479754
## 17 Rima 1.2587930
## 18 Rua_a 1.8735826
## 19 Tahi_a 2.4392563
## 20 Tahi_c 1.6405735
## 21 Toru 1.3542783
## 22 Whaa 0.9798064
mean(a$timelaps)## [1] 2.490691
# what's the maximum time between locations?
a <- aggregate(timelaps~segmentid, raw_tracks, max)
a## segmentid timelaps
## 1 212499 22.39028
## 2 212500_a 21.43056
## 3 215258_a 70.58194
## 4 215258_b 36.63500
## 5 215259_a 49.24583
## 6 215261_a 22.82778
## 7 215261_b 37.23389
## 8 215262_a 34.74778
## 9 215262_b 21.46528
## 10 215262_d 51.34000
## 11 215263 22.48583
## 12 46633 21.35639
## 13 46635 42.47000
## 14 46950 18.19167
## 15 46955 10.59694
## 16 Bill Wiremu 23.38111
## 17 Rima 73.75361
## 18 Rua_a 20.35194
## 19 Tahi_a 22.30056
## 20 Tahi_c 15.07583
## 21 Toru 61.70472
## 22 Whaa 10.33222
save(raw_tracks, file = "./Outputs/raw_tracks.RData")g <- gbase +
geom_sf(data = raw_tracks, aes(col = segmentid, geometry = geometry), size = 0.8)+
theme(legend.position = "right")
ggsave(g, file = "./Outputs/filtered_tracks_NZ_segments.png", dpi = 300, width = 8, height = 8)
gFigure 1.3: Filtered satellite tracks of right whales tagged in the Auckland Islands in 2020 and 2021 and South Georgia in 2020. Any segment of less than 15 positions separated by interruptions of more than 6 days was removed
Prepare data in a specific format for foieGras package
raw_tracks$time <- raw_tracks$Date
data_ssm <- raw_tracks %>%
dplyr::select(segmentid, time, Quality, Longitude, Latitude, Error.Semi.major.axis, Error.Semi.minor.axis, Error.Ellipse.orientation) %>%
dplyr::rename(id = segmentid, date = time, lc = Quality, lon = Longitude, lat = Latitude, smaj = Error.Semi.major.axis, smin = Error.Semi.minor.axis, eor = Error.Ellipse.orientation)
data_ssm$geometry <- NULL
str(data_ssm)## 'data.frame': 32329 obs. of 8 variables:
## $ id : chr "212499" "212499" "212499" "212499" ...
## $ date: POSIXct, format: "2021-07-21 09:26:18" "2021-07-21 09:29:49" ...
## $ lc : chr "A" "B" "A" "B" ...
## $ lon : num 166 166 166 166 166 ...
## $ lat : num -50.5 -50.5 -50.5 -50.5 -50.5 ...
## $ smaj: int NA 1110 1751 6119 1359 1769 2128 4656 4081 5304 ...
## $ smin: int NA 222 43 117 802 1332 892 1109 1121 758 ...
## $ eor : int NA 61 106 172 166 93 117 93 98 107 ...
save(data_ssm, file = "./Outputs/data_ssm.RData")Fit spate-space model (based on codes from Reisinger et al. 2021)
nz_fit_ssm <- tibble(time_step = c(3, 6)) %>% # create a tibble with the two time step conditions in one column
group_by(time_step) %>% #
nest() %>%
mutate(data = list(data_ssm, data_ssm), # store the same dataframe in the two rows
fit = map2(time_step, data, function(x, y){
fit_ssm(y,
spdf = F, # we already run the filter so no need to do again
pf = F,
model = "rw",
min.dt = 48,
time.step = x,
control = ssm_control(verbose = 0),
map = list(psi = factor(NA)))}))Assess movement persistence along tracks
nz_fit_ssm <- nz_fit_ssm %>%
mutate(fmp = map(fit, function(a){
fit_mpm(a, what = "predicted", model = "mpm", control = mpm_control(verbose = 0))
}))
save(nz_fit_ssm, file = "./Outputs/nz_fit_ssm.RData")load("./Outputs/nz_fit_ssm.RData")
pdf(file = "./Outputs/nz_fit_ssm_interpolation.pdf", width = 12, height = 12)
nz_fit_ssm %>%
pull(fmp) %>%
llply(., function(x) {plot(x, pages = 1, ncol= 5, pal = "Zissou1", rev = T)})## [[1]]
##
## [[2]]
dev.off()## png
## 2
nz_fit_ssm %>%
pull(fmp) %>%
llply(., function(x) {plot(x, pages = 1, ncol= 5, pal = "Zissou1", rev = T)})## [[1]]
Figure 2.1: Time series of movement persistence modeled along right whale tracks from Auckland Islands. The move persistence parameter gt ranges continuously from 0 (little persistence, indicative of area-restricted movements) to 1 (high persistence, indicative of directed movements).
##
## [[2]]
Figure 2.2: Time series of movement persistence modeled along right whale tracks from Auckland Islands. The move persistence parameter gt ranges continuously from 0 (little persistence, indicative of area-restricted movements) to 1 (high persistence, indicative of directed movements).
Convert predictions to a dataframe and extract ARS
nz_fit_ssm_df <- nz_fit_ssm %>%
unnest(fmp) %>%
dplyr::select(-c(fit, data)) %>%
mutate(pred = map(mpm, function(x){
cbind(x$data[c("date", "x", "y")], x$fitted["g"])
})) %>%
dplyr::select(c(time_step, id, pred)) %>%
unnest(pred) %>%
st_as_sf(coords = c("x", "y"), crs = 4326) %>%
st_transform(proj_antarctic)
# distribution of g values
summary(nz_fit_ssm_df$g)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.003713 0.418424 0.640095 0.587120 0.771512 0.976121
# select the 50% quantile as a threshold for ARS. it matches pretty well the cut-off used by Zerbini et al. with the beta parameter of the SSM
nz_fit_ssm_df <- nz_fit_ssm_df %>%
mutate(mode = case_when(g < quantile(g, 0.5) ~ "ARS",
g >= quantile(g, 0.5) ~ "transit"))fun_plot_ars <- function(ids) {
for(i in ids){
d <- nz_fit_ssm_df %>% filter(id == i)
g1 <- ggplot(d, aes(col = g)) +
geom_sf() +
facet_wrap(~time_step) +
scale_color_viridis_c(name = "Move persistence")
g2 <- ggplot(d, aes(col = mode)) +
geom_sf() +
geom_sf(data = d[d$mode == "ARS", ], aes(col = mode)) + # plot again on top
facet_wrap(~time_step)
g <- (g1 / g2) &
plot_annotation(title = i)
print(g)
}
}
pdf(file = "./Outputs/nz_fit_ssm_interpolation_tracks.pdf", width = 10, height = 8)
fun_plot_ars(unique(nz_fit_ssm_df$id))
dev.off()## png
## 2
fun_plot_ars(unique(nz_fit_ssm_df$id))Figure 4.1: Movement persistence modeled along right whale tracks from Auckland Islands. The move persistence parameter gt ranges continuously from 0 (little persistence, indicative of area-restricted movements) to 1 (high persistence, indicative of directed movements).
Figure 4.2: Movement persistence modeled along right whale tracks from Auckland Islands. The move persistence parameter gt ranges continuously from 0 (little persistence, indicative of area-restricted movements) to 1 (high persistence, indicative of directed movements).
Figure 4.3: Movement persistence modeled along right whale tracks from Auckland Islands. The move persistence parameter gt ranges continuously from 0 (little persistence, indicative of area-restricted movements) to 1 (high persistence, indicative of directed movements).
Figure 4.4: Movement persistence modeled along right whale tracks from Auckland Islands. The move persistence parameter gt ranges continuously from 0 (little persistence, indicative of area-restricted movements) to 1 (high persistence, indicative of directed movements).
Figure 4.5: Movement persistence modeled along right whale tracks from Auckland Islands. The move persistence parameter gt ranges continuously from 0 (little persistence, indicative of area-restricted movements) to 1 (high persistence, indicative of directed movements).
Figure 4.6: Movement persistence modeled along right whale tracks from Auckland Islands. The move persistence parameter gt ranges continuously from 0 (little persistence, indicative of area-restricted movements) to 1 (high persistence, indicative of directed movements).
Figure 4.7: Movement persistence modeled along right whale tracks from Auckland Islands. The move persistence parameter gt ranges continuously from 0 (little persistence, indicative of area-restricted movements) to 1 (high persistence, indicative of directed movements).
Figure 4.8: Movement persistence modeled along right whale tracks from Auckland Islands. The move persistence parameter gt ranges continuously from 0 (little persistence, indicative of area-restricted movements) to 1 (high persistence, indicative of directed movements).
Figure 4.9: Movement persistence modeled along right whale tracks from Auckland Islands. The move persistence parameter gt ranges continuously from 0 (little persistence, indicative of area-restricted movements) to 1 (high persistence, indicative of directed movements).
Figure 4.10: Movement persistence modeled along right whale tracks from Auckland Islands. The move persistence parameter gt ranges continuously from 0 (little persistence, indicative of area-restricted movements) to 1 (high persistence, indicative of directed movements).
Figure 4.11: Movement persistence modeled along right whale tracks from Auckland Islands. The move persistence parameter gt ranges continuously from 0 (little persistence, indicative of area-restricted movements) to 1 (high persistence, indicative of directed movements).
Figure 4.12: Movement persistence modeled along right whale tracks from Auckland Islands. The move persistence parameter gt ranges continuously from 0 (little persistence, indicative of area-restricted movements) to 1 (high persistence, indicative of directed movements).
Figure 4.13: Movement persistence modeled along right whale tracks from Auckland Islands. The move persistence parameter gt ranges continuously from 0 (little persistence, indicative of area-restricted movements) to 1 (high persistence, indicative of directed movements).
Figure 4.14: Movement persistence modeled along right whale tracks from Auckland Islands. The move persistence parameter gt ranges continuously from 0 (little persistence, indicative of area-restricted movements) to 1 (high persistence, indicative of directed movements).
Figure 4.15: Movement persistence modeled along right whale tracks from Auckland Islands. The move persistence parameter gt ranges continuously from 0 (little persistence, indicative of area-restricted movements) to 1 (high persistence, indicative of directed movements).
Figure 4.16: Movement persistence modeled along right whale tracks from Auckland Islands. The move persistence parameter gt ranges continuously from 0 (little persistence, indicative of area-restricted movements) to 1 (high persistence, indicative of directed movements).
Figure 4.17: Movement persistence modeled along right whale tracks from Auckland Islands. The move persistence parameter gt ranges continuously from 0 (little persistence, indicative of area-restricted movements) to 1 (high persistence, indicative of directed movements).
Figure 4.18: Movement persistence modeled along right whale tracks from Auckland Islands. The move persistence parameter gt ranges continuously from 0 (little persistence, indicative of area-restricted movements) to 1 (high persistence, indicative of directed movements).
Figure 4.19: Movement persistence modeled along right whale tracks from Auckland Islands. The move persistence parameter gt ranges continuously from 0 (little persistence, indicative of area-restricted movements) to 1 (high persistence, indicative of directed movements).
Figure 4.20: Movement persistence modeled along right whale tracks from Auckland Islands. The move persistence parameter gt ranges continuously from 0 (little persistence, indicative of area-restricted movements) to 1 (high persistence, indicative of directed movements).
Figure 4.21: Movement persistence modeled along right whale tracks from Auckland Islands. The move persistence parameter gt ranges continuously from 0 (little persistence, indicative of area-restricted movements) to 1 (high persistence, indicative of directed movements).
Figure 4.22: Movement persistence modeled along right whale tracks from Auckland Islands. The move persistence parameter gt ranges continuously from 0 (little persistence, indicative of area-restricted movements) to 1 (high persistence, indicative of directed movements).
g <- gbase +
geom_sf(data = nz_fit_ssm_df, aes(col = g, geometry = geometry), size = 0.8) +
theme(legend.position = "right") +
scale_color_viridis_c(name = "Move persistence") +
facet_wrap(~time_step)
ggsave(g, file = "./Outputs/nz_fit_ssm_tracks_all.png", dpi = 300, width = 16, height = 8)
gFigure 4.23: Interpolated satellite tracks at 1 position/3hours of right whales tagged in the Auckland Islands in 2020 and 2021.
g <- gbase +
geom_sf(data = nz_fit_ssm_df, aes(col = mode, geometry = geometry), size = 0.8) +
geom_sf(data = nz_fit_ssm_df[nz_fit_ssm_df$mode == "ARS",], aes(col = mode, geometry = geometry), size = 0.8) +
theme(legend.position = "right") +
facet_wrap(~time_step)
ggsave(g, file = "./Outputs/nz_fit_ssm_tracks_ars.png", dpi = 300, width = 16, height = 8)
gFigure 4.24: Movement mode of right whales tagged in the Auckland Islands in 2020 and 2021.